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ABSTRACT 

We present results from simulations of rotating magnetized turbulent convection in spherical wedge geometry 
representing parts of the latitudinal and longitudinal extents of a star. Here we consider a set of runs for 
which the density stratification is varied, keeping the Reynolds and Coriolis numbers at similar values. In the 
case of weak stratification we find quasi-steady solutions for moderate rotation and oscillatory dynamos with 
poleward migration of activity belts for more rapid rotation. For stronger stratification a similar transition as a 
function of the Coriolis number is found, but with an equatorward migrating branch near the equator. We test 
the domain size dependence of our results for a rapidly rotating run with equatorward migration by varying 
the longitudinal extent of our wedge. The energy of the axisymmetric mean magnetic field decreases as the 
domain size increases and we find that anm = l mode is excited for a full 2ir 0-extent, reminiscent of the field 
configurations deduced from observations of rapidly rotating late-type stars. 

Subject headings: Magnetohydrodynamics - convection - turbulence - Sun: dynamo, rotation, activity 



1. INTRODUCTION 

The large-scale magnetic field of the Sun, manifested by 
the 1 1 year sunspot cycle, is generally believed to be gener- 
ated within or ju st below the turbulent convection zone (e.g. 
lOssend riiver 2003, and references therein). The latter con- 
cept is based on the idea that the strong shear in the tachocline 
just beneath the convection zone amplifies the toroidal mag- 
netic field which the n becomes buo yantly unstable and erupts 
to the surface (e.g. Parker 1955b). This process has been 
adopted in many mean-field model s of the solar cycle in the 
form of a non-local a-effect (e.g. Kitchati nov & Ole mskov 
120121) in conjunction with a reduced turbulent diffusivity 
within the convection zone, and a single cell anti-clockwise 
meridional circulation which acts as a conveyor belt for 
the magnetic field. These so-call ed flux transport models 
(e.g. iDikpati & Charbonneaulll999l) are now widely used to 
study the solar cycle and even to predict its f uture course 
(Dikp ati & Gilmadl27x^rChoudhuri et alll2007h . 

The flux transport paradigm is, however, facing several the- 
oretical challenges: 10 5 ga uss magnetic fields are ex pected 
to reside in the tachocline (iD'Silva & C houdhurij ll993l) . but 
such fields are difficult to explain with dynamo theory 
dGuerrero & Kapvla1l201 ll) a nd may have bec ome unstable at 
much lower field strengths (lArlt et al.l 120051) . Furthermore, 
the low value of turbulen t diffusion within the convection 
zone seems unlikely (e.g. IKapyla et alJl2009h . In addition, 
the meridional circulation pattern of the Sun can be more 
complicated than that required in the flux-transport models 
dHalhawavll20TltlMiesch et alj|2012l) . These dif ficulties have 
led to a revival of the distributed dynamo (e.g. iBrandenburgl 
12005b iPipinl 120121) in which magnetic fields are generated 
throughout the convection zone due to turbulent effects (e.g. 
Krause & RadlellHM IKapyla et alJ[200rj : iPipin & Seehafed 



2009). 



Early studies of self-consistent three-dimensional magne- 
tohydrodynamic (MHD) simulations of convection in spher- 
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ical coordinates produced oscil latory large-scale dynamos 
(Gilmar] [l983t iGlatzmaier 1985), but the dynamo wave was 
found to propagate towards the poles rather than the equator 
as in the Sun. More recent anelastic large-eddy simulations 
(LES) with rotation rates somewhat higher than that o f the 
Sun have p roduced non-oscillatory (|Brown et all 120101) and 
oscillatory dBrown et al.ll201 lllNelson et al.ll2013l) large-scale 
magnetic fields, depending essentially on the rotation rate and 
the vigor of the turbulence. However, similar models with the 
solar rotation rate have e ither failed to pr oduce an appreciable 
large-scale component dBrun et all 120041) or, more recently, 
oscillatory solutions with almost no lat i tudinal propagatio n 
of the activity belts dGhizaru et alj|2010t iRacine et al.ll201 ll) . 
These simulations covered a full spherical shell and used re- 
alistic values for solar luminosity and rotation rate, neces- 
sitatin g the use of an elastic solvers and spheri cal harmonics 
fe.g. l Brun et al.ll2004l) or implicit methods (e.g. lGhizaru et al.l 
2010). Here we exploit an alternative approach by modeling 
fully compressible conv ection in wedge geometry (see also 
Robin son & Cha n 2001) with a finite-difference method. We 
omit the polar regions and cover usually only a part of the 
longitudinal extent, e.g. 90° instead of the full 360°. At the 
cost of omitting connecting flows across the poles and in- 
troducing artificial boundaries there, the gain is that higher 
spatial resolution can be achieved. Furthermore, retaining 
the sound waves can be beneficial when considering pos- 
sible h elio- or asteroseismi c application s. Recent hydrody- 
namic (IKapyla et alJl2011allbh and MHD dKapyla et alll 201(1) 
studies with this method have shown that this approach pro- 
duces results that are in accordance with fully spherical mod- 
els. Moreover, the first turbulent dynamo solution with solar- 
like migration properties of the magnet ic field was re cently 
obtained using this type of setup (Kap yla et alJl2012bT) . Ex- 
tended setups that include a coronal layer as a more real- 
istic upper radial boundary have bee n successful in produc - 
ing dynamo-driven coronal ejections (Warn ecke et al.ll2012l) . 
As we show in a companion paper ( Warnec ke et al.ll2013l) . a 
solar-like differential rotation pattern might be another conse- 
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quence of including an outer coronal layer. 

Here we concentrate on exploring further the recent dis- 
covery ofeguatowardjriigration in spherical wedge simula- 
tions (iKapyla et alJl2012bh . In particular, we examine a set 
of runs for which the rotational influence on the fluid, mea- 
sured by the Coriolis number, which is also called the inverse 
Rossby number, is kept approximately constant while the den- 
sity stratification of the simulations is gradually increased. 

2. THE MODEL 

Our model is the same as that in Kapv la et alj d2012al) . We 
model a wedge in spherical polar coordinates, where (r, 9, <fr) 
denote radius, colatitude, and longitude. The radial, latitudi- 
nal, and longitudinal extents of the wedge are r < r < R, 
9q < 6 < it — 9o, and < <fi < <fio, respectively, where R is 
the radius of the star and ro = 0.7 R denotes the position of 
the bottom of the convection zone. Here we take 9q = tt/12 
and in most of our models we use <fio = tt/2, so we cover 
a quarter of the azimuthal extent between ±75° latitude. We 
solve the compressible hydromagnetics equations 1 , 

dA 

= u x B r]^oJ, (1) 



dt 



Du 
~Dt 



D\i\p 
Dt 
1 



= -V • it, 



(2) 



g - 2£l xu + -(J x B -Vp + V ■ 2vpS) , (3) 



T- 



Dt 



P 



rad 



pSGS^ 



L M) r ] J 2 ]+2vS 2 , (4) 



where A is the magnetic vector potential, u is the velocity, 
B = V x A is the magnetic field, J = /i ( 7 1 V x B is the 
current density, r\ is the magnetic diffusivity, /io is the vacuum 
permeability, D/Dt = d/dt + u ■ V is the advective time 
derivative, v is the kinematic viscosity, p is the density, 



rirad 



KVT and F 



SGS 



XsgspTVs (5) 



are the radiative and subgrid scale (SGS) heat fluxes, where K 
is the radiative heat conductivity and xsgs turbulent heat con- 
ductivity, which represents the unresol ved convective trans - 
port of heat and was referred to as xt in Kap vla et alj d2012al) . 
s is the specific entropy, T is the temperature, and p is the 
pressure. The fluid obeys the ideal gas law with p = (7— l)pe, 
where 7 = cp/cy = 5/3 is the ratio of specific heats at con- 
stant pressure and volume, respectively, and e = cyT is the 
internal energy. 

The gravitational acceleration is given by g = —GMrjr 2 , 
where G is the gravitational constant, M is the mass of 
the star, and r is the unit vector in the radial direction. 
In simulations, the maximum possible Rayleigh number is 
much smaller than in real stars, which implies a larger 
Mach number. The rotation vector CIq is given by Oo = 
J7o(cos 9, — s'md, 0). However, to have realistic Coriolis 
numbers, the angular velocity in the Coriolis force is in- 
creased correspondingly, but that in the centrifugal force is 
omitted, as it woul d otherwise be unrealistically large (cf. 
Kapvla et al. 2 01 lbl) . The rate of strain tensor S is given by 



(6) 



1 Note that in Equation (4) of Kapyla' et al. 1 2012b) the Ohmic heating term 
fior/J 2 and a factor p in the viscous dissipation term 2^S 2 were omitted, but 
they were actually included in the calculations. 



where the semi colons denote covariant differentiation 
dMitra et al.ll2009l) . 

2. 1. Initial and boundary conditions 

The initial state is isentropic and the hydrostatic tempera- 
ture gradient is given by 

dT _ GM/r 2 

dr Cy(7 — l)(n + l)' 

where n = 1.5 is the poly tropic index. We fix the value of 
dT/dr on the lower boundary. The density profile follows 
from hydrostatic equilibrium. The heat conduction profile is 
chosen so that radiative diffusion is responsible for supplying 
the energy flux in the system, with K decreasin g more than 
two or ders of magnitude from bottom to top (IKapyla et al.1 
1201 lab . A weak random small-scale seed magnetic field is 
taken as initial condition (see below). 

The radial and latitudinal boundaries are assumed to be im- 
penetrable and stress free, i.e., 

du e u g <9u u 
u r = 0, ^- = — , = — (r = r , R), (8) 

or r 

8u r <9U0 

For the magnetic field we assume perfect conductors at the 
lower radial and latitudinal boundaries, and radial field at the 
outer radial boundary. In terms of the magnetic vector poten- 
tial these translate to 

° X * 1 - (r = r ), (10) 



dr 

■ cot ( 



(e = e ,7T-e ). (9) 



dr 

A r = 0, 



= Ag = = 
dAg Ag 



dr 



dA^ 
dr 



A d 



%,7r-e ). 



(r = R), (11) 



(12) 



We use small-scale low amplitude Gaussian noise as initial 
condition for the magnetic field. On the latitudinal bound- 
aries we assume that the thermodynamic quantities have van- 
ishing first derivatives, thus suppressing heat fluxes through 
the boundaries. 
On the upper boundary we apply a black body condition 



dr 



dr 



(13) 



where a is the Stefan-Boltzmann constant. We use a modi- 
fied value for a that takes into account that our Reynolds and 
Rayleigh numbers are much smaller than in reality, so K and 
therefore the energy flux through the domain are much larger 
than in the Sun. 

2.2. Dimensionless parameters 

To facilitate with other work using different normalizations, 
we present our results in this paper by normalizing with phys- 
ically meaningful quantities. We note, however, that in the 
code we used non-dimensional quantities by choosing 



R = GM = po 



cp 



Mo 



1. 



(14) 



where po is the initial density at r — ro. The units of length, 
time, velocity, density, entropy, and magnetic field are there- 
fore 



[x] = R, [t] = ^/R 3 /GM, [u] = y/GM/R, 
p = Po , [s] = cp, [B] = y/poVLoGM/R. (15) 
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Our simulations are defined by the energy flux imposed at the 
bottom boundary, Fb = — (KdT/dr)\ r=ro , the temperature 
at the top boundary, T\ = T(r = R), as well as the values of 
Oq> v, an d Xsgs = XSGs(?"m = 0.85 R). Furthermore, the 
radial profiles of K and xsgs need to be specified. The corre- 
sponding nondimensional input parameters are the luminosity 
parameter 

C= Po (GAf)3/2i?i/2' (16) 
the normalized pressure scale height at the surface, 

(7 - l)cvli 



1.00 



^ GM/R ' 
the Taylor number 

Ta = (2Q a R 2 /iy) 2 , 
the fluid and magnetic Prandtl numbers 



Pr = 



Pm = -, 



XSGS 

and the non-dimensional viscosity 

v 



VGMR 



(17) 



(18) 



(19) 



(20) 



Instead of £, we often quote the resulting density contrast, 
r p = p(r )/p(R). 

Other useful diagnostic parameters are the fluid and mag- 
netic Reynolds numbers 



Re = 



Mid 



vkf ' 



Pan = 



rjk{ ' 



and the Coriolis number 

CO: 



2Q 

ZZrms&f 



(21) 



(22) 



where u lms = J (3/2)(u 2 . + Ug) r g^,t is the rms velocity and 
the subscripts indicate averaging over 9, <fi, and a time interval 
during which the turbulence is statistically stationary. Note 
that for u rms we omit the contribution from the azimuthal ve- 
locity, because its value is dominated b y effects from the dif- 
ferential rotation (Kapyl iiet al.ll201 lbl) . The Taylor number 
can also be written as Ta = Co 2 Re 2 (fcf-R) 4 , with kfR w 21. 
Due to the fact that the initial stratification is isentropic, we 
quote the turbulent Rayleigh number Rat from the thermally 
relaxed state of the run, 



Ra t 



GM(Ar) 4 



1 d(s) g< p t 
cp dr 



(23) 



where k{ — 2n/Ar is an estimate of the wavenumber of the 
largest eddies, and Ar = R — tq = 0.3R is the thickness of 
the layer. We also quote the value of fc w = uj lms /u lma , where 
u> = V x u and u) lma , is the volume averaged rms value 
of u>. The magnetic field is expressed in equipartition field 

strengths, B cq (r) — (popu 2 )^^., where all three components 
of u are included. We define mean quantities as averages over 
the ^-coordinate and denote them by overbars. However, as 
we will see, there can also be significant power in low-order 
spherical harmonic modes with azimuthal order m = 1 and 2, 
which will be discussed at the end of the paper. 
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FIG. 1. — Initial (solid lines) and saturated (dashed) radial profiles of tem- 
perature T, density p, and pressure p, normalized by their respective values 
at the bottom of the domain from Run CI (indicated by subscript zero). The 
inset shows the specific entropy s/cp from the same ran. 

The simulations were performed with the PENCIL Code 2 , 
which uses a high-order finite difference method for solving 
the compressible equations of magnetohydrodynamics. 

3. RESULTS 

We perform runs for four values of £, corresponding to ini- 
tial density stratifications T p = 2, 5, 30, and 100. These runs 
are referred to as series A-D. In series E we use T — 30 and 
vary <j>o with all other parameters being kept the same as in 
Run CI. For each series we consider different values of Co 
and Re. The hydrodynamic progenitors of the Runs B 1 , CI, 
andDl correspond to Runs A4, B4, and C4 from Kapy la et alj 
( 1201 lab . The rest of the simula tions were run from the initial 
conditions described in Section lXTl 

Earlier studies applying fully spherical simulations have 
shown that organized large-scale magnetic field s appear pro- 
vided the rotation of the star is rapid enough (Bro wn et al.l 
2010) and that at even higher rotation r ates, cyclic soluti ons 
with poleward migration are obtained dBrown et al.112.01 ll) . A 
simila r transition has been obse rved in spheri c al wed ge mod- 
els of lKapvla etaT] d20lob and lKapvla etal] <2012aD . How- 
ever, in the former case the oscillatory mode showed poleward 
migration of the activity belts whereas in the latter an equa- 
torwar d branch appears near the equator. Furthermore, in the 
runs of Kapvla et alj d2012al) the dynamo mode changes from 
one showing a high frequency cycle with poleward migration 
near the equator to another mode with lower frequency and 
equatorward migration when the magnetic field becomes dy- 
namically important. 

There are sever al d ifferences be t ween the models of 
iKapvla et al] ( 120101) and lKapvla etal] ( 1201 2al) : the amount of 
density stratification (a density contrast of 3 in comparison 
to 30), efficiency of convective energy transport (20 per cent 
versus close to 100 per cent in the majority of the domain), 
and the top boundary condition for entropy (constant temper- 
ature versus black body radiation). Here we concentrate on 
studying the influence of th e density strat i fication on models 
similar to those presented in Kap vla et alj d2012al) . 

3.1. Thermal boundary effects and energy balance 

In IKapvla et alj ( 1201 1 ah we started to apply the black 
body boundary condition, Equation ( fT3l , that has previ- 
ously been used in mean-field models with thermodynamics 

2 http://code.google.eom/p/pencil-code/ 
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TABLE 1 
Summary of the runs. 



Run 


grid 


Pr 


Pm 


Ta 


c 






V 






C 




Ra t 


Re 


Rm 


Co 


A 1 


1ZO X ZOO X IzCS 


1.5 


1.0 


1 nlO 
10 


0.29 


2 


1.7 


10~ 


5 


3.8 


10" 


-5 


8.3 


10 5 


26 


26 


8.6 


A2 


128 X 256 X 128 


1.5 


1.0 


1.8 ■ 10 


0.29 


2 


1.7 


10~ 


5 


3.8 


10~ 


5 


1.1 


10 5 


24 


24 


12.8 


Bl 


128 x 256 x 128 


2.5 


1.0 


6.4 • 10 a 


0.09 


5 


2.9 


10~ 


b 


3.8 


io- 


!> 


1.1 


10 b 


22 


22 


8.1 


B2 


128 x 256 x 128 


2.5 


1.0 


1.4 ■ 10 10 


0.09 


5 


2.9 


10~ 


5 


3.8 


io- 


5 


1.1 


10 6 


20 


20 


13.7 


CI 


128 x 256 x 128 


2.5 


1.0 


1.4 • 10 1U 


0.02 


30 


2.9 


10~ 


b 


3.8 


io- 


b 


2.1 


10 b 


35 


35 


7.8 


C2 


128 x 256 x 128 


2.5 


1.0 


4.0 ■ 10 10 


0.02 


30 


2.9 


io- 


5 


3.8 


io- 


5 


2.7 


10 6 


31 


31 


14.8 


Dl 


128 x 256 x 128 


7.5 


3.0 


1.6 • 10 a 


0.008 


100 


4.7 


io- 


6 


6.3 


io- 


<; 


1.2 


10 b 


11 


34 


8.0 


D2 


256 x 512 x 256 


4.0 


2.0 


10 10 


0.008 


100 


2.5 


io- 


5 


6.3 


io- 


6 


2.4 


10 6 


25 


50 


9.1 


El 


128 x 256 x 64 


2.5 


1.0 


1.4 ■ 10 1U 


0.02 


30 


2.9 


10~ 


5 


3.8 


10" 


5 


2.1 


10 b 


34 


34 


7.9 


E2 


128 x 256 x 128 


2.5 


1.0 


1.4 ■ 10 10 


0.02 


30 


2.9 


io- 


5 


3.8 


io- 


5 


2.1 


10 6 


35 


35 


7.8 


E3 


128 x 256 x 256 


2.5 


1.0 


1.4 ■ 10 10 


0.02 


30 


2.9 


io- 


5 


3.8 


io- 


5 


2.4 


10 6 


35 


35 


7.9 


E4 


128 x 256 x 256 


2.5 


1.0 


1.4 ■ 10 10 


0.02 


30 


3.5 


io- 


5 


3.8 


io- 


5 


2.2 


10 6 


28 


28 


8.1 



NOTE. — The second to ninth columns show quantities which are input parameters to the models whereas the quantities in the last four columns are results 
of the simulations computed from the saturated state. Here we use 0o — tt/2 in Sets A-D. In Set E we use (ftp — 7r/4 (R un El), (ftp — 7r/2(E2),^o — ^ 
(E3), and <pp — 2tt (E4). Runs CI and E2 are the same model, and is also the same as Run B4m of Kiipyla et al. i 2012a). 



TABLE 2 

Summary of diagnostic variables. 



Run 


A 


-Emcr/i?kin 


-Erot/-Ekin 


-Emag / Ekin 


Bpol/ ^mag 


Etor 1 Em&g 




A (9) 




Al 


0.084 


0.000 


0.580 


0.418 


0.045 


0.396 


0.013 


0.089 


62 


A2 


0.095 


0.000 


0.490 


0.553 


0.068 


0.338 


0.009 


0.050 


62 


Bl 


0.028 


0.000 


0.705 


0.345 


0.038 


0.487 


0.034 


0.142 


68 


B2 


0.098 


0.000 


0.757 


0.222 


0.056 


0.427 


0.023 


0.072 


72 


CI 


0.006 


0.001 


0.440 


0.346 


0.138 


0.203 


0.047 


0.068 


93 


C2 


0.105 


0.001 


0.326 


0.706 


0.198 


0.238 


0.016 


0.030 


94 


DI 


0.003 


0.002 


0.222 


0.472 


0.166 


0.135 


0.011 


-0.000 


89 


D2 


0.003 


0.000 


0.617 


0.222 


0.133 


0.190 


0.045 


0.058 


116 


El 


0.007 


0.001 


0.478 


0.393 


0.133 


0.328 


0.048 


0.069 


92 


E2 


0.006 


0.001 


0.440 


0.346 


0.138 


0.203 


0.047 


0.068 


93 


E3 


0.005 


0.001 


0.375 


0.380 


0.120 


0.172 


0.037 


0.055 


92 


E4 


0.024 


0.001 


0.363 


0.470 


0.065 


0.114 


0.035 


0.050 


91 



NOTE. — Here A — A/ (u ims fcf ) is the normalized growth rate of the magnetic field. Skin — \ (p u2 ) is the volume averaged 
kinetic energy. fi mer — ^(p{u^. -f u'$)) and E TO t — ^{p^e) denote the volume averaged energies of the azimuthally averaged 
meridional circulation and differential rotation. Analogously E ma , s — \{B 2 ) is the total volume averaged magnetic energy while 
.Epoi — i ((B 2 , + B 2 ))) and Et or — \ are tne energies in the axisymmetric part of the poloidal and toroidal magnetic fields. 



F am v = cp((pu r )'T'), (25) 
F kin = -±((pu r yu' 2 ), (26) 

F v i sc = -2v{(pUi)' Sir) , (27) 

^turb = -XSGs(p)(T)^, (28) 

Fpoyn = (E e B^ - EfB )/no, (29) 

where E = tj/j-qJ u x B, the primes denote fluctuations, 
and angle brackets averaging over 9, <f>, and a time interval 
over which the turbulence is statistically stationary. The ra- 
diative flux carries energy into the convection zone and drops 
steeply as a function of radius so that it contributes only a few 
per cent in the middle of the convection zone. The resolved 
convection is responsible for transporting the energy through 
the majority of the layer, whereas the unresolved turbulent 
transport carries energy through the outer surface. The other 
contributions are much smaller. The flux of kinetic energy is 
also ver y small in the rapid ro tation regime considered here 
(see also[Augustso n et al.ll2012h . 

3.2. Dynamo excitation and large-scale magnetic fields 

The azimuthally averaged toroidal magnetic fields from all 
runs listed in Table [T] are shown in Figures |3]-|6] The full time 



((Brandenburg et al. 1992). Instead of using the physical value 
for the Stefan-Boltzmann constant, we estimate the value of 
a so that the flux at the upper boundary is approximately that 
needed to transport the total luminosity of the star through 
the surface. However, the final thermally relaxed state of the 
simulation can significantly deviate from the initial state. In 
combination with the nonlinearity of Equation ( fT3l l, the fi- 
nal stratification is usually somewhat different from the initial 
one; see Figure [TJ for an illustrative example from Run CI. 
The final density stratification in this case is around 22, down 
from 30 in the initial state. 

The main advantage of the black body condition is that it 
allows the temperature at the surface more freedom than in 
our previous models where a constant temperature is imposed 
(Kapvlaet al. 2010, 2011b). Furthermore, as the temperature 
is allowed to vary at the surface, this can be used as a diag- 
nostic for possibl e irr adiance variations. These are discussed 
further in Section [3~6l 

Considering the energy balance, we show the energy fluxes 
for Run E4 in Figure |2| We find that the simulation is ther- 
mally relaxed and that the total luminosity is close to the input 
luminosity, i.e. i t ot — Lq ~ 0. The fluxes are defined as: 

- d(T) 

F iad = -K-^, (24) 
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FIG. 2. — Luminosity of the energy fluxes from Run E4: radiative conduc- 
tion (thin solid line), enthalpy (dashed), kinetic energy (dot-dashed), viscous 
(long dashed), and unresolved subgrid scale (dotted) fluxes. The thick solid 
line is the sum of all contributions. The two dashed red lines indicate the zero 
and unity line. 

evolution from the introduction of the seed magnetic field 
to the final saturated state is shown for each run. Although 
the magnitude of the seed field in terms of the equipartition 
strength is different in each set, already a visual inspection 
suggests that the growth rate of the dynamo is greater in the 
cases with low stratification. The measured average growth 
rates over the kinematic stage, 

X=(dlnB ims /dt) u (30) 

support this conjecture; see the second column of Table [2] 
Comparing Runs Al, Bl, CI, and D2 with roughly compara- 
ble Reynolds and Coriolis number shows that the normalized 
growth rate decreases monotonically from 0.084 in Run Al 
to just 0.003 in Run D2. Another striking feature is that A 
increases by a factor of nearly 20 from Run CI to C2 where 
the only difference between the runs is that the latter has a 
roughly two times higher Coriolis number. It turns out that in 
all of the cases (Runs Al, A2, Bl, B2, C2, and E4) with the 
highest growth rates, a poleward migrating dynamo mode at 
low latitudes is excited first. In some of the runs this mode 
is later overcome by another one that can be quasi-stationary 
(Runs Al and Bl) or oscillatory with equatorward migration 
and a much longer cycle period (Runs C2 and E4). 

In Figure|3]we show the azimuthally averaged toroidal mag- 
netic field near the surface of the computational domain 
(r = 0.98 R) for two runs (Al and A2; see Table [T} with 
T p = 2 (Set A). We find that in Run Al with Co 8.7 
the mean magnetic field is initially oscillatory with poleward 
propagation of the activity belts. At tu lma k{ « 400 the 
dynamo mode changes to a quasi-steady configuration. In 
Run A2 a poleward mode persists throughout the simula- 
tion, although the oscillation period is irregular and significant 
hemispherical asym metry exists. This b ehavior is similar to 
the runs presented in Ka pyla et al.l (12010) which are compara- 
ble in terms of stratification, Coriolis and Reynolds numbers. 

In Set B with T p = 5 the situation is similar: in Run Bl 
with Co w 8.1 there is a poleward mode with a short cycle 
period near the equator which is visible from early times, see 
Figure H] However, after around tu lms k{ = 1200 there is a 
dominating non-oscillatory mode which is especially clear at 
high latitudes. There are still hints of the poleward mode near 
the equator. In Run B2 with Co w 13.7, however, the pole- 
ward mode prevails also at late times. As in Run A2, the cy- 
cles show significant variability and hemispheric asymmetry. 
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FIG. 3. — near the surface of the star at r = 0.98 R as a function of 
latitude (= 90° - 9) for Runs Al (top) and A2 (bottom). The white dotted 
line denotes the equator 90° —9 = 0. 

The runs in Sets A and B a lso show signs of non-axisymme tric 
'nests' of convection (cf. Busse 2002; Brown et al. 2008) in 
the hydrodynamical and kinematic stages. Once the magnetic 
field is strong enough these modes either vanish or they are 
significantly damped. 

Increasing the stratification further to T p = 30 (Set C) the 
dynamo solutions at lower rot ation rates, Co < 5, are still 
quasi-steady; see Figure 2 of iKapvla et alj d2012al) . How- 
ever, a watershed regarding the oscillatory modes at higher 
Co seems to have been reached so that the irregular poleward 
migration seen in Sets A and B is replaced by more regular 
equatorward patterns. In Run CI with Co w 8.7 the poleward 
migration near the equator is secondary to the equatorward 
mode — even in the kinematic stage; see Figure [5] The pole- 
ward mode near the equator is more prominent in the early 
stages of Run C2 with Co ps 14.7, but also there it is subdom- 
inant at late times. 

For r p = 100 (Set D) the general picture is similar to that 
in Set C. Quasi-steady configurations at lower rotation change 
into equatorward migrating solutions at sufficiently high Co. 
We find this transition to occur between Co = 5 and 8, simi- 
larly as in Set C. For Set D the equatorward mode is visible for 
both runs; see Figure @] In Run Dl no poleward migration at 
low latitudes is seen in the kinematic stage. Also the poleward 
migrating branch at high latitudes is missing in the non-linear 
stage. Both of these features are present in Run D2. 

To connect our results with observations of magnetically 
active stars we compute the ratio of the dynamo cycle fre- 
quency and the rotation rate, cj cvc /Oo. Plotting this ra- 
tio as a function of the Coriolis number for stars exhibit- 
ing chromospheric activity has shown that stars tend to 
group along inactive and active branches (Brandenburg et al. 
1998), and for higher Coriolis num bers along a super-active 
branch dSaar & Brandenburg 1999). Six of our simulations 
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FIG. 5. — Same as Figure[3]but for Runs CI (top panel) and C2 (bottom). 
Note the difference in cycle frequency between early times when the fre- 
quency is similar to that of Run B2 (Figure[4j and late times. 



(Runs A2, B2, CI, C2, Dl, and D2), excluding the runs in 
Set E which are very similar to each other and Run CI in this 
respect, show cycles and can be thus used in this analysis. For 
Runs CI, C2, Dl, and D2 the cycle frequency was measured 
from a few latitudes on both hemispheres and for all available 



cycles individually. The final result is an average over the in- 
dividually computed values. In two of our models (Runs A2 
and B2) the cycle periods vary irregularly. There the cycle 
frequency is computed from the highest peak of a temporal 
Fourier transformation of the time series for at low lat- 
itudes near the surface. The results are shown in the upper 
panel of Figure [7J 

The association with real branches in Figure [7J is clearly 
premature, because we cannot be sure that there are no stars 
connecting the group of Runs CI, Dl, D2 with that of A2 
and B2 through a single line with a larger slope. Neverthe- 
less, this plot allows us to see that, while the separation in the 
ratio w cyc / £!o is slightly less for the two groups of runs, com- 
pared with active and inactive stars, their relative ordering in 
the value of Co is actually the other way around. One would 
therefore not have referred to Runs A2 and B2 as inactive just 
because their w cyc /Oo ratio agrees with that of inactive stars. 
In fact, their E mag / 'Skin ratios in Tableware typically larger 
than for Runs CI, Dl, D2. 

As is clear from the inset of the upper panel of Figure |7J 
there is no clear relation between Co and E mag / E^, which 
is clearly different from stars for which there is a clear re- 
lation between Co (referred to as inverse Rossby number in 
tha t context) and stellar activ ity (a measure of -Bmag/^kin); 
see iBranden burg et al. ( 1998) for details and references. Fur- 
thermore, there is also no clear indication of two branches in 
the graph of oj cyc /flo versus E mag / E^; see the lower panel 
of Figure [TJ Instead, there might just be one group in it, pos- 
sibly with a positive correlation, i. e., uj rvc /^o might increas e 
with Smag/^kin- As discussed by Brandenburg et al. ( 1998), 
a positive slope would not be easy to explain in the framework 
of standard mean-field dynamo theory, where the frequency 
ratio is usually a decreasing function of normalized rotation 
rate and activity parameter. 

In conclusion, we can say that the quantity ui cyc /flo is an 
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FIG. 7. — Ratio of the dynamo cycle and rotation frequencies from six 
runs which show cyclic activity as functions of rotation (upper panel) and 
the ratio of magnetic to kinetic energy (lower panel). The dotted and dashed 
lines in the up per are given by CiCo"' , where the ai correspond to those in 
Brandenburg et al. ( 1998), and Cj are used as fit parameters. The inset in the 
top panel show the ratio Emag/^Hn as a function of Co. 





FIG. 8. — Mean rotation profile f2/f2o (contours) and meridional circulation 
u m = (u r ,ug, 0) (aiTows) from Sets A, B, C, and D. 

important and robust property of any cyclic dynamo model 
and its dependence on other properties of the model should 
therefore be a useful characteristics that can be compared with 
other models and ultimately with actual stars. Here we have 
just make a first attempt in classifying model results in this 
sense. 



3.3. Differential rotation and meridional circulation 

Non-uniform rotation of the convection zone of the Sun is 
an important ingredient in maintaining the large-scale mag- 
netic field. Furthermore, the sign of the radial gradient of 
the mean angular velocity plays a crucial role in deciding 
whether the dynamo wave propagat es towards the pole o r the 
equator in mean-field models (e.g. Parker 1955a! [1987I) . In 
the following we use the local angular velocity defined as 
£1 = ul/r sin 9 + f2rj. Azimuthally averaged rotation profiles 
from the runs in Sets A to D are shown in Figure [8] The ro- 
tation profiles of Runs El, E3, and E4 are very similar to that 
of Run CI. We quantify the radial and latitudinal differential 
rotation by 



(r) ^eq ^bot 

^ = n ' 

(6) ^cq ^pole 



n 



(31) 
(32) 



cq 



where fl cq = £l(R, 7r/2) and f2bot = Q( r o, 7r /2) are the an- 
gular velocities at the top and bottom at the equator, respec- 
tively, and n polc = ffi(R, O ) + fi(-R, tt - $ )]. 
Mean-field and LES models of stellar rotation are most of- 



ten performed without including magnetic fields (e.g.j lempel 
12001 iKitchatinov & Rudigeill2 005: Mi esch et alj|2006h . This 
can be misleading because magnetic fields affect the tur- 
bule nce giving rise to Reynolds stress and a heat flux 
(e.g. IKitchatinov etai] 11994 IKapvla et alj [2004). Further- 
more, the large-scale flows are directly influenced by the 
Lorentz force when the magnetic field is strong enough (e.g. 

Malk us~& Proctor! 1 19751) . A de crease of * ha s also been 
observed in LES models (e.g. iBrun et al.l 12004). Compar- 
ing the latitudinal differential rotation in Run Bl with that 
of otherwise identical hydrodynamic Run A4 of Kapy la et al.1 

d201 lal) . we find that Afo' decreases only slightly from 0.15 

(r) 

to 0. 14. For Aj 2 ; the change is more dramatic from 0.079 to 
0.034. The fraction of kinetic energy contained in the differ- 
ential rotation drops from 0.91 to 0.71. A similar decrease 
is observed in Run CI in comparison to its hydrodynamical 

parent Run B4 of IKapvla et alj d2011al) with A^ changing 

from 0.08 to 0,07, A^ from 0.066 to 0.047, and E mt /E kin 
dropping from 0.58 to 0.44. 

We see a rapidly spinning equator with a positive radial gra- 
dient of il in all cases; see Figure [8] The latitudinal vari- 
ation of angular velocity is however not always monotonic 
and there can be local minima at mid-latitudes, as is seen for 
example i n Run CI. Similar features have been s een before 
(see, e.g., lMiesch et alj|20 00; Kapvla et al. 2011b) and might 
be related to the lack of small-scale turbulence. Especially at 
larger stratification one would expect smaller-scale turbulent 
structures to emerge, but this requires sufficient resolution and 
thus large enough Reynolds numbers, which is not currently 
possible. 

The amount of latitudinal differential rotation is clearly 



less than in the Sun where A 



(0) 



.2 bet ween the equa- 



tor and latitude 60° (e.g. ISchou et alj[l998l) . Furthermore, 

A n generally decreases within each set of runs as Co in- 
creases with the exception of Runs Dl and D2. However, 
in Run Dl the lower Reynolds number possibly contributes 
to the weak differential rotation in comparison to Run D2 
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FIG. 9. — Dynamo number C a from Sets A, B, C, and D. 

with comparable Co. The rotation profiles appear to be dom- 
inated by the Taylor-Proudman balance except at very low 
latitudes wher e the baroclinic term is significant. In a com- 
panion paper (Warn ecke et al.ll2013h . we show that an outer 
coronal layer seems to favor a solar-like rotation, which shows 
even radially outward pointing contours of constant rotation. 
Such 'spoke-like' rotation profiles have so far been obtained 
only in mean-field models invo lving anisotropic heat transport 
(e.g. Brandenburg et al. 1992) or a subadiabatic tachocline 
(Rempel 2005), and in LES models where a lat itudinal en- 
tropy gradient is enforced at the lower boundary (Mies ch et alj 
2006). 

The meridional circulation is weak in all cases and typically 
shows multiple cells in radius. The circulations are mostly 
concentrated in the equatorial regions outside t he inner tan- 
gent cy linder. This is similar to earlier results by Kapv la et al.l 
(2012a) where the meridional circulation pattern was shown 
in larger magnification. In addition, as we will show below, 
the strength of meridional circulation relative to the turbu- 
lent magnetic diffusivity is rather low, which is another reason 
why it cannot play an important role in our models. 

3.4. Estimates of dynamo numbers 

We estimate the dynamo numbers related to the a-effect, 
radial differential rotation, and the meridional circulation by 

dTl/dr(AR) 3 



FIG. 10. — Dynamo number Cq from Sets A, B, C, and D. We omit regions 
closer to 2.5° from the latitudinal boundaries. 



aAR 

Vto 



Cc 



a 



iAR 



Vto 



Vto 



(33) 

where d^l/dr is the r- and ^-dependent radial gradient of £1 
and a is a proxy of the a-effect 



a = -|r(a; u-j- b/p), 



(34) 



with t = «MLT-ffp /u rms (r, 8) being the convective turnover 
time and ckmlt the mixing length parameter. We use »mlt = 
5/3 in this work. We estimate the turbulent diffusivity 
by vto = ^u rms (r,6)a MLT H P . Furthermore, m™. = 

y/v!i + u'g is the rms value of the meridional circulation. 

The results for the dynamo numbers are shown in Fig- 
ures [914TT1 Generally, the values of C a are fairly large, and 
those of Co surprisingly small, suggesting that the dynamos 
might mainly be of or type. In the following, however, we 




FIG. 1 1 . — Dynamo number Cjj from Sets A, B, C, and D. 



focus on relative changes between different runs. It turns out 
that there is a weak tendency for C Q to increase as a function 
of r p and Co. In Set B, however, C a decreases by a third 
from Run Bl to B2. The spatial distribution of C a becomes 
more concentrated near the boundaries as T p increases. 

We find that the effect of the differential rotation is strongest 
near the equator in all cases. In Sets A and B the absolute 
values of Co are clearly smaller that those of C Q . This is 
surprising given the fact that the energy of the toroidal mean 
field is greater than that of the poloidal field by a significant 
factor which would be expected if differential rotation domi- 
nates over the a-effect in maintaining the field. In Runs CI, 
C2, and Dl C a and Co have comparable magnitudes whereas 
in Run D2 the maximum of C a is roughly twice that of Cq. 
However, in these cases the toroidal and poloidal field ener- 
gies are roughly comparable. 

We find that Cu is always small in comparison to both C a 
and Co- Figure [TT] also shows the concentration of coherent 
meridional circulation cells in the equatorial regions with a 
multi-cell structure. 



FIG. 12. — Radial velocity u r (top row) and toroidal magnetic field (bottom) near the surface of the star r = 0.98 R in Mollweide projection from a Runs El 
(left), E2, E3, and E4 (right). 



3.5. Ejfect of domain size 

We recently reported equatorward migration of activity 
belts in a spherical wedge simulation (Kap yla et al.ll2012al) . 
There we gave results from simulations with a ^-extent of 
tt/2. However, at large values of the Coriolis number, 
the a effect becomes sufficiently anisotropic and differen- 
tial rotation weak so that nonaxisym metric solutions become 
possible; see Moss & Brandenburg d 19951) for corresponding 
mean-field models with dominant m = 1 modes in the limit 
of rapid rotation. To allow for such modes, we now choose a 
0-extent of up to 2tt. In the present case, we find that for Co ps 
7.8 it is possible that non-axisymmetric dynamo modes of low 
azimuthal order (to = 1 or 2) can b e excited. Th i s woul d 
not be possible in the simulations of Kapvla et aj] d2012al) . 
The same applies to non-axis ymmetric modes excited in hy- 
drodynamic convection (e .g. Busse 2002; B rown et~aT1l2008t 
iKapyla et aT]|201 lbtlAugustson et alJl2012h ~ 

We test the robustness of this result by performing runs with 
A0 = 7r/4, 7r/2, 7T, and 2ir with otherwise identical parame- 
ters. We find that the same dynamo mode producing equator- 
ward migration is ultimately excited in all of these runs. The 
only qualitatively different run is that with c/)q = 2tt where the 
poleward mode near the equator grows much faster than in the 
other cases. However, after tu rms k{ s» 1500 the equatorward 
mode takes over similarly as in the runs with a smaller <pQ. 

The velocity field shows no marked evidence of non- 
axisymmetry but there are indications of m = 1 structures 
in the instantaneous magnetic field; see Figure Q~2] This is 
also reflected by the fraction of the axisymmetric component 
of the total magnetic energy, see the 5th and 6th columns of 
Table |2] We find that the energy of the mean toroidal field 
decreases monotonically when 0o is increased so that there 
is a factor of three difference between the extreme cases of 



Runs El and E4. The axisymmetric component still exhibits 
an oscillatory mode with equatorward migration in all runs in 
Set E. 

We compute power spectra of the toroidal component of 
the magnetic field from the Run E4 over three 10° latitude 
strips from each hemisphere, centered around latitudes ±25°, 
±45° and ±65°. The results for the three lowest modes are 
shown in Figure Qj] We find that at low (±25°) and high 
(±65°) latitudes the axisymmetric (to = 0) mode begins 
to dominate after around 1000 turnover times and shows a 
cyclic pattern consistent with that seen in the time-latitude di- 
agram of the azimuthally averaged field. After tu lms k{ ~ 
1600, however, the m = 1 mode becomes stronger in the 
southern hemisphere. This coincides with the growth of 
the to = 1 mode at mid-latitudes where it dominates al- 
ready earlier in both hemispheres. This is in rough agree- 
ment with some observational results of rapid rotators, which 
show the most prominent non-axisymmetric t e mperature (e.g. 
Hackmanetal. 2001; Korhonen et al. 2007; Lindborg et al. 
201 ll) and magnetic structures (Kochuk hov et al.l l2013) at the 
latitudinal range around 60-80°, while the equatorial and po- 
lar regions are more axisymmetric; some temperature inver- 
sions even show almost completely axisymmetric distribu- 
tions in the p olar regions and ri ngs of azimuthal field at low 
latitudes (e.g. lDonati et al.l2003l) . The strength of the axisym- 
metric versus the non-axisymmetric component has also been 
reported to vary over ti me with a time scale of a few years 
dKochukhov et alj|20l3h . 

3.6. Irradiance variations 

The black body boundary condition allows the temperature 
vary at the surface of the star and thus enables the study of ir- 
radiance variations due to the magnetic cycle. In Figure[T4lwe 
show for Run CI the difference between the azimuthally av- 
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FIG. 13. — Energies of the m = (black lines), 1 (red), and 2 (blue) 
modes of the toroidal magnetic field as functions of time near the surface of 
the star (r = 0.98 R) in Run E4. The data is averaged over 10° latitude 
strips centered at latitudes 90° — 8 = ±25 (top panels), ±45 (middle), and 
±65 (bottom) degrees and normalized by the total energy within each strip. 

eraged temperature and its temporal average for the saturated 
state of the dynamo. We find a clear signal which is strongest 
for the poleward branch at high latitudes but also visible for 
the equatorward branch near the equator. The peak values at 
high latitudes are up to 20 per cent of the surface tempera- 
ture. This is relativel y large compared with ea rlier work using 
mean-field models (Brandenburg et al. 1992), which showed 
remarkably little relative variation of the order of 10~ 3 in the 
bulk of the convection zone and even less at the surface. This 
difference is probably related to the importance of latitudinal 
variations that were also present in the mean-field model of 
Brandenbu rg et al. l (11992) and referred to as thermal shadows 
(lParkerlll987h ~ 

4. CONCLUSIONS 

We have studied the effects of density stratification on the 
dynamo solutions found in simulations of rotating turbulent 
convection in spherical wedge geometry. For the four values 
of r p , which is the ratio of the densities at the bottom and 
at the surface of the convection zone. In addition, we vary 
the rotation rate for each value of T p . For all stratifications 
we find quasi-steady large-scale dynamos for lower rotation 
and oscillatory solutions when rotation is rapid enough. The 
transition from quasi-steady to oscillatory modes seems to oc- 
cur at a lower Co for higher stratification. Furthermore, for 
low values of T p the oscillatory solutions show only poleward 
propagation of the activity belts whereas in higher T p an equa- 
torward branch appears at low latitudes. 

The equatorward branch was first noticed by Kapyla et al. 
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FIG. 14. — Temperature fluctuation at the surface from Run CI. The values 
are saturated at 7 per cent of the mean, whereas the maxima at high latitudes 
are around 20 per cent. The white dotted line denotes the equator. 



d2012al) using a wedge with = 90° longitude extent. Here 
we test the robustness of this result by varying Acf> from 45° 
to full 360°. We find a very similar pattern of the axisymmet- 
ric component of the field in all cases. However, the energy 
of the axisymmetric magnetic field decreases with increasing 
Atfi. In the simulation with the full 2ir ^-extent we observe 
an m = 1 mode which is visible even by visual inspection 
of the simulation data (see Figure [T2l . Such field configu- 
rations have been observed in rapidly rotating late-type stars 
(see e.g. Kochu khov et alj2013l) and our simulation is the first 
time such features have been seen in direct numerical simula- 
tions. We are currently investigating the rapid rotation regime 
with more targeted runs which will be reported in a separate 
publication (Cole et al. 2013, in preparation). 

The ratio between cycle to rotation frequency, Lu cyc /fl , is 
argued to be an important non-dimensional output parameter 
of a cycle dynamo. For the Sun and other relatively inac- 
tive stars, this ratio is around 0.01, while for the more ac- 
tive stars it is around 0.002. For our models we find values 
in the range 0.002-0.01, but for most of the runs it is around 
0.004. Although it is premature to make detailed comparisons 
with other stars and even the Sun, it is important to emphasize 
that kinematic mean-field dynamos produce the correct cycle 
frequency only for values of the turbulent magnetic diffusiv- 
ity that are at least 10 times smaller tha n what is suggested 
by standard estimates (IChoudhurilll990l) . In our case, these 
longer cycle periods (or smaller cycle frequencies) are a re suit 
of nonlinearity and are only obtained in the saturated regime 
of the dynamo. The detailed reason for this is unclear, but it 
has been speculated that it is connecte d with a slow magnetic 
helicity evolution (Brandenburg 2005). Equally unclear is the 
reason for equatorward migration, which, as we have seen, is 
a consequence of nonlinearity, too. It will therefore be impor- 
tant to provide an accurate determination of all the relevant 
turbulent transport coefficients. 
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